#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2018 crane <crane@crane-pc>
#
# Distributed under terms of the MIT license.

"""

"""

from matrix_tools import *
from gaussian_eliminate import Gaussian
from matrix_class import Matrix

g = Gaussian()
M = Matrix()

class Proj:
    def __init__(self, onto_mat):
        ''' 默认投影到onto_mat的列空间 '''
        self.proj_onto = self._col_basis(onto_mat)
        self.proj_onto_t = transpose_matrix(self.proj_onto)

        self.t_multiply_proj = multiply(self.proj_onto_t, self.proj_onto)

        self.proj_mat = self.project_mat()


    def _col_basis(self, onto_mat):
        d = M.four_space_basis(onto_mat)
        a = d['col']
        return transpose_matrix(a)

    def project_mat(self):
        ''' 求投影矩阵
            A* (AT*A)^-1 * AT

            if A's column is not independent, extract independent column.
            get rid of redundant column
        '''
        a = self.proj_onto
        aT = self.proj_onto_t
        i = g.I(self.t_multiply_proj)

        return multiply(a, i, aT)

    def project(self, v):
        return multiply(self.proj_mat, v)

    def x(self, v):
        ''' onto_mat * x = proj_v

                 AT * (b - A * x) = 0
            ===> AT * A * x = AT * b
            ===> x = (AT * A)^-1 * AT * b
        '''
        a = self.proj_onto
        aT = self.proj_onto_t
        i = g.I(self.t_multiply_proj)
        return multiply(i, aT, v)


def test():
    # 故意构造列相关矩阵(计算投影矩阵时, 会先获取column basis)
    # 这里的[1, 1, 1] 是常量C, 或者说是截距.
    onto = [
        [1, 1],
        [2, 1],
        [3, 1],
        [4, 1],
    ]
    pj = Proj(onto)
    # pj = Proj([
    #     [1],
    #     [2],
    #     [3],
    # ])
    # pj = Proj([
    #     [2],
    #     [4],
    #     [6],
    # ])

    p = pj.proj_mat
    show_matrix(p, 'proj mat')

    # 投影矩阵平方还是本身
    p2 = matrix_power(p, 2)
    # show_matrix(p2, 'proj square')

    # 投影矩阵对称
    # assert is_symmetric([[1,2], [2,3 ], [3,4]])
    assert is_symmetric(p)
    assert is_symmetric(p2)

    # =================== proj something  =====================
    v = [
        [1],
        [3],
        [8],
        [9],
    ]

    project_v = pj.project(v)
    show(project_v, 'pj vector')

    trans, rref = M.rref(onto)
    show(trans, 'trans')
    show(rref, 'rref')

    rref_for_v = multiply(trans, project_v)
    show(rref_for_v, 'rref_for_v')

    # =================== proj x =====================
    x = pj.x(v)
    show(x, 'x')

def test2():
    ''' 满秩矩阵的投影矩阵是identity单位矩阵 '''
    pj = Proj([
        [1, 1, 2],
        [2, 3, 4],
        [3, 2, 7],
    ])

    p = pj.proj_mat
    show_matrix(p, 'proj mat')

    # 投影矩阵平方还是本身
    p2 = matrix_power(p, 2)
    show_matrix(p2, 'proj square')

    # 投影矩阵对称
    # assert is_symmetric([[1,2], [2,3 ], [3,4]])
    assert is_symmetric(p)
    assert is_symmetric(p2)
    assert is_identity(p2)

    # =================== proj something  =====================
    v = [
        [1],
        [3],
        [5],
    ]

    project_v = pj.project(v)
    show(project_v, 'pj vector')


def main():
    test()
    # test2()

if __name__ == "__main__":
    main()
